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By using a non-local model, fluid simulations can capture kinetic effects in the parallel electron heat-flux better 
than is possible using flux limiters in the usual diffusive models. Non-local and diffusive models are compared using 
a test case representative of an ELM crash in the JET SOL, simulated in one dimension. The non-local model 
shows substantially enhanced electron temperature gradients, which cannot be achieved using a flux limiter. The 
performance of the implementation, in the BOUTH — h framework, is also analysed to demonstrate its suitability for 
application in three-dimensional simulations of turbulent transport in the SOL. 



The divertor target heat-flux will be an important limit 
on the performance of future magnetic confinement fusion 
devices, from ITER onwards[l, 2]. Modelling of trans- 
port in the scrape-off layer (SOL) is therefore a sub- 
ject of urgent interest. The conditions in the SOL — 
large amplitude fluctuations, open magnetic field lines — 
prevent the application of approximations that are used 
to make (gyro-)kinetic simulations of turbulence feasible 
in the core plasma; kinetic simulations are limited to 
one-dimensional models (e.g. [3]). Nevertheless, three- 
dimensional features — turbulence and 'blob '[4] motion — 
are critical parts of the picture, especially in determin- 
ing the width of the scrape-off layer [5, 6]. These are 
studied using fluid simulations which are either two- 
dimensional (e.g. ESEL[7, 8], RI/FI-S0L[9], SOLT[10] and 
G-ESEL[11]) or use simple (Braginskii) models for the par- 
allel dynamics (e.g. heat diffusion in GBS[12] or Spitzer 
conductivity in B0UT++[13, 14]). However, such simple 
treatments of the parallel transport may not give the full 
picture. As a consequence of an edge localized mode 
(ELM) or the motion of a blob through the SOL there 
will be large changes in the overall density and temper- 
ature on a field line. When such transient events occur 
kinetic simulations show that the diffusive description is 
not a good model of the parallel transport: the correc- 
tion for 'kinetic effects' that can be applied by using flux 
limiters breaks down since the appropriate limiter, as de- 
termined by comparison with kinetic simulations, varies 
strongly in time [15]. 

To enable a self-consistent treatment of turbulence with 
'kinctically-corrected' parallel transport, new methods for 
including kinetic effects in fluid models are required. This 
paper describes the application of a non-local calculation 
of the electron heat-flux [16], which has not previously been 
applied to the SOL. The parallel electron heat-flux is de- 
rived by using a very high order truncation (with typically 
100 to 1600 moments retained) to solve a one-dimensional 
reduced drift kinetic equation in terms of a set of integrals. 
This model has a couple of particularly attractive features. 
Firstly, being derived directly from the kinetic equation, 
it has no ad-hoc parameters. Secondly, a good deal of 
the computational complexity involved is found in solving 
for the eigensystem that is used to decouple the moment 



equations. This need be done only once for a given trun- 
cation, and so a substantial piece of the 'kinetic' part of 
the calculation is removed from the simulations. Alternat- 
ives, which have been used in laser-plasma physics[17, 18], 
employ integral kernels which would be numerically very 
challenging when, as is the case here, the density cannot 
be assumed to be constant. 

The model has been implemented in the B0UT++ 
plasma fluid simulation framework[19] in order to expedite 
its future inclusion in three-dimensional simulations. Here 
we evaluate its performance through one-dimensional sim- 
ulations of parallel transport, with parameters chosen to 
approximate an ELM in JET[20]. 

In Section 1 we outline the model and detail some addi- 
tions needed for its use in the SOL; in Section 2 the results 
of simulations of the test case are presented, comparing the 
non-local model to the Braginskii model with and without 
flux limiters; in Section 3 the numerical performance of 
the implementation is discussed, which is particularly im- 
portant for future extensions to three-dimensional fluid 
models; finally in Section 4 we conclude and discuss fu- 
ture developments. 

1 Theory 

1.1 Non-local heat-flux 

The thermal speed of the electrons in SOL plasmas is very 
much higher than that of the ions. Hence in calculations 
of the electron dynamics, we may consider the evolution 
of the background profiles to be slow. Solving a static 
electron kinetic equation for the heat-flux will therefore 
give a good approximation to its true value, which can 
then be used in the evolution equations for the background 
fields. 

The approach taken in [16] starts from the one- 
dimensional kinetic equation for the non-Maxwellian part 
5f e of the electron distribution function whose Maxwellian 
part is fe°^'- 
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where (•) denotes the gyroaverage, vu is the component of 
the velocity parallel to the magnetic field and C (•, •) is the 
linearized Fokker-Planck collision operator. We expand 

in fluid moments on a basis P ifc (v) = P l (v)L k 2 (v 2 ), 

P'(v) being tensor harmonic polynomials and L^ +2 \v 2 ) 
associated Laguerre polynomials (see [21] for details), and 
truncate to L angular harmonics and K Laguerre orders, 
i.e. < I < L, < k < K; here we will take L = 20, K = 
20 and so have 400 moments in total. After truncating 
it is convenient to represent the pairs (I, k) with a single 
index, A,B,..., running over Lx K possible values. Then 
(1) reduces to a set of one-dimensional, first order ODEs 
for the non-Maxwellian fluid moments 
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where z is the dimensionless length defined by || = ^- 
(with Ac the electron-electron collision length) ; n A are the 
parallel fluid moments of f e ; ^ A B is the matrix coming 
from v\ 1 ; C B is the collision matrix; and 
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is the drive term from the gradient of the Maxwellian 
part, fe . These equations are valid only for the non- 
Maxwellian moments (i.e. A, B = (0,0), (0,1) and (1,0) 
are excluded): the Maxwellian moments (the density, tem- 
perature and fluid velocity) are background fields which 
give the drive term g A . 

The term due to the gradient of the thermal velocity, 
which was neglected in [16], has been retained; as we will 
see below it is unimportant for the core plasma, but it 
becomes relevant in the SOL. The equations (2) can be 
decoupled using the eigenvectors of the matrix operator 
(C* _1 *) A B to give 
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where n A and g A are the components of n A and g A on the 
eigenvector basis and Qa) are the corresponding eigen- 
values (note that as defined here Qa) are the reciprocals 
of the eigenvalues used in [16] 1 ). The solutions of these 

1 The reason for this discrepancy is that in [16] the eigenvalues 
of (\E'~ 1 C) j4 s are found. Due to the block-off-diagonal structure 
of 'J', when the three Maxwellian moments are removed the upper- 
left-most blocks are no longer square, with the result that ^ is not 
invertible. The authors of [16] resolved this by changing the trun- 
cation to keep three extra moments (0, K), (0, K + 1) and (1, K) so 
that all the blocks of remain square. C, though, is block-diagonal 
and so is invertible with or without the extra moments. The ei- 
genvalues and eigenvectors of (C" 11 !/) g can therefore be found in 
either case, though without the extra moments there will be (ex- 
actly one) zero eigenvalue. The approach taken here is thus more 
flexible. In this paper we do not add the three extra moments, and 
do have a zero eigenvalue. The effect is to change the response to 
fluctuations whose wavelength is short enough (compared to the col- 
lision length) that the model does not converge. Without a zero 
eigenvalue the heat-flux due to such fluctuations is anomalously low, 
while with one it is anomalously high. The latter case may be better 
behaved if, for instance, numerical noise happened to introduce such 
short-wavelength fluctuations. 



equations can be expressed as 
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where fi A (z a ) are the values at the boundary z (where 
zq = Z- < z for ((A) < and zq = z + > z for ((A) > 0). 
The factors result from the term. In the core, 



flux surfaces are closed and so fluctuations in T are small 
compared to T, thus t(zo) ~ Tjz^j ~ ^ anc ^ tnese factors 
can be neglected (as they are in [16]). In the SOL this is 
not the case as flux surfaces are open and the change in 
temperature along a field line may be large. 

Finally, after transforming back to the original vector 
basis we use the (1,1) moment to calculate the parallel 
heat-flux, which is needed to close the fluid equations: 
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where W A B — is the matrix formed from the eigen- 

vectors W( B ) (whose eigenvalues are C(s))- 

So to calculate g e || we must (a) compute integrals over 
z (as many as the number of moments retained) and (b) 
set the boundary values n A (zo) (which we will do so as to 
impose the boundary condition on the heat-flux, Section 
1.3). 

1.2 Fluid Equations 

For our simulations we use the Id fluid equations, assum- 
ing quasineutrality and ambipolarity, 
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where = (■^ + ^ / ||V||), S n is the particle source, S B is 
the energy source (which is taken to be equal for electrons 
and ions) and r e i is the electron-ion collision time. The 
equations are closed by specifying q e ||, and 71411: 

• <7 e 1 1 is given either by the non-local model described 
above, by the Braginskii diffusion operator q c i. e = 
— 3.16 " e ^ eT " V||T e , or by the flux- limited diffusion op- 
erator 
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Figure 1: Target heat-flux for simulations using non-local electron heat-flux (left) and unlimited diffusive electron 
heat-flux (right). The electron component is the red dashed line, the ion component is the blue dotted line and the 
total is the black solid line. N.B. this is the target heat-flux (at the material surface) and not the sheath-edge heat-flux 
that is used as the boundary condition on the fluid model, i.e the electron heat-flux has been shifted down and the ion 
heat-flux shifted up by the sheath potential. The heat-flux has not, however, been adjusted for the incidence angle at 
the target, so if the angle were, for instance, 6° the true heat-flux at the material surface would be smaller by a factor 
of - 0.1. 



where the free-streaming heat-flux is \qt 8 e | = nT e VTe 
and a e is a parameter to be set. 

• gj| | is calculated using a flux-limited diffusion operator 

9i\\ = (— + ^—) (12) 

where the Braginskii diffusion operator is q^ = 
-3.9 " T ^ T " V||T t , the free-streaming heat-flux is 
k/s,i| = nTiVTi and on is a parameter to be set. 

• 7Tj|| is given by a limited form of the Braginskii vis- 
cosity 

"" = (4^) s J0rfj (13) 

where ir = — | x .96nTiy/2TuV \\ Vj| and b is a para- 
meter to be set. 

1.3 Boundary Conditions 

In the SOL the boundary conditions are a critically im- 
portant part of the dynamics. Here we wish to impose the 
sheath-edge fluid boundary conditions[22] 

V se = ±c s = ±f T ^ fJ ± (14) 
V m 

q se . e = (2T e + \e<p s \) nc s (15) 

where ecf> s = 0.5T e In (l + is the sheath po- 

tential and 7 = 3, representing an approximately collision- 
less sheath. 

(14) and (15) under-determine the boundary values of 
the non-Maxwellian moments h A . Since our aim here is 
to develop computationally efficient methods to improve 
fluid turbulence models, we simply choose the boundary 



moments so as to impose (15) on the heat-flux and leave 
the remaining moments unchanged, by adding a contribu- 
tion proportional to (W ) ^ jy 

There is an additional complication, that in a hot 
enough plasma (which is found in the simulations de- 
scribed below) the collision length may be long enough 
that the transients originating at one boundary do not de- 
cay to zero by the time they reach the other. This means 
that the boundary conditions cannot be set locally: one 
requires information from both boundaries to set them 
consistently. 

2 Simulations 

We conducted simulations to compare the non-local elec- 
tron heat-flux model with the usual diffusive model. The 
test case chosen is that used by [20], which is represent- 
ative of the JET SOL. An ELM crash is simulated by 
introducing transient sources of heat and particles with 
duration 200/is whose spatial distributions are cosines 
about the centre of the simulation domain. The total 
ELM energy is 0.4MJ, which is distributed evenly between 
the ions and electrons, both at the pedestal temperature 
Tped = 1.5keV. The source region of the SOL is taken to 
have a width of 10cm, a radius of 3m, a poloidal extent of 
2.6m and a field line pitch of 6°, giving a volume of 4.9m 3 
and a parallel extent of 2.6m/ sin (6°) ~ 25m. Explicitly 
the source terms are 




where y is the co-ordinate parallel to the magnetic field 
with y = at the centre of the domain, L s = 25m is 
the length of the source, a n « 9.1 x 10 23 m _3 s _1 and 
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Figure 2: Electron temperature profiles at various times during and after a 200/is ELM using the non-local heat-flux 
model (left) and a diffusive heat-flux model with a flux limiter of 0.4 (right). 
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Figure 3: Comparison of the electron temperature profiles from the non-local (solid, blue) and diffusive models with 
a flux limiter of 0.4 (dashed, red) and without flux limiter (dotted, green): on the left at t = 0/j,s (before the ELM) 
and on the right at t — 50/iS (during the ELM). 



<je ~ 3.3 x 10 8 Jm _3 s _1 . The total length of the domain is 
80m. The initial profiles are found by allowing the simula- 
tion to relax to a steady state with sources that have sim- 
ilar form but smaller amplitude (<j n « 7.4 x 10 22 m _3 s _1 
and oe ~ 4.5 x 10 6 Jm _3 s _1 ). The ion heat-flux limiter 
used was ctfj =0.1 and the ion viscosity limiter was b = 0.5, 
as indicated by inter-ELM PIC simulations [15]. There 
are 256 uniformly distributed spatial grid points and a 
staggered grid is used, with fluxes evaluated at the cell 
edges. 

It has been noted before[3, 20] that the target heat-flux 
is not strongly affected by kinetic effects, the most im- 
portant factor being the ions whose transport is mostly 
convective and therefore well described by fluid models. 
It is not surprising then that the target heat-flux is sim- 
ilar for the non-local and the diffusive heat-flux models, 
as shown in Figure 1. One might have thought that us- 
ing a model which includes kinetic effects on the electron 
heat-flux would show the rapid response due to fast elec- 
trons that is seen in PIC simulations[20]. This is not 
the case here because both source terms and boundary 
conditions are given in the fluid picture and do not in- 



clude kinetic corrections. In order to give a description 
more fully in agreement with the PIC results some modi- 
fication of the fluid boundary conditions would also be 
necessary: in particular the sheath electron heat trans- 
mission coefficient shows very strong variation through an 
ELM[15]. One proposal to achieve this is to use sheath 
coefficients corresponding to the response to a collisionless 
Maxwellian wavepacket[23], but this relies on knowledge 
of the source to derive an explicit time dependence. For 
three-dimensional simulations where the source on a par- 
ticular field line is not known in advance, a method which 
depends only on the fields being evolved is needed. The 
non-local calculation here depends only on the fluid vari- 
ables, but derives from them extra information about the 
non-Maxwellian part of the electron distribution function. 
It would be interesting to investigate whether this can 
be put to use in the construction of kinetically-corrected 
boundary conditions. The conclusion of this is that for 
one-dimensional studies the non-local electron heat-flux 
does not by itself offer much improvement with respect 
to modelling of machine-relevant parameters, which are 
dominated by the ion dynamics. 
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The value of the non-local model lies instead in as- 
pects which will be relevant when coupled in to three- 
dimensional simulations. As shown in Figures 2 and 3 the 
non-local model has a strong effect on the shape of the 
electron temperature profiles. The qualitative difference 
between this non-local model and a flux limiter approach 
is clear. Before the ELM, a flux limiter with a e — 0.4 
is in fairly good agreement with the non-local model, but 
during the ELM, when the electron temperature is much 
higher and the collision length is much longer, the same 
flux limiter cannot sustain temperature gradients of a sim- 
ilar order of magnitude as the non-local model. To do so a 
much smaller flux limiter would be needed, for example at 
t = 50/is a limiter with a e ~ 0.025 would be required (see 
Figure 4). A flux limiter cannot be set which is appropri- 
ate for the whole range of conditions. This is in agreement 
with PIC simulations of the SOL[f 5], which show that to 
capture the kinetic effects the parameter of a flux limiter 
model would need to vary both in space and time. In 
contrast the non-local heat-flux can cover a much broader 
range of conditions. A PIC simulation of the same sort of 
ELM (0.4MJ, 200/zs) as we have simulated here shows[15, 
Fig. 7] that the value of the electron heat-flux limiter that 
would match the kinetic simulation becomes much smaller 
than the pre-ELM value after 50/zs, showing that in the 
kinetic simulation a much steeper temperature gradient is 
needed to drive the heat-flux from the source than would 
be the case with a flux limiter set to the pre-ELM value. 
Figure 4 shows a similar plot from our non-local heat-flux 
simulation: the spatial average, (a e ), of the flux limiter 
parameter (here un-normalized, unlike [15, Fig. 7]) that 
would be needed at each point for the diffusion operator 
(11) to match the non-local heat-flux. In both cases the 
minimum value reached by (a e ) is similar. [15, Fig. 6] 
shows that the kinetic correction to the electron sheath 
heat transmission coefficient, j e , passes through 1 (i.e. no 
correction) at about 60/is, which happens to be the time 
when (a e ) reaches its minimum value. This is why our 
non-local fluid simulation, which does not include any kin- 
etic corrections to the fluid boundary conditions, is able to 
find the same value. After 60/is, the fluid heat transmis- 
sion coefficient is larger than the kinetic one, and so the 
electrons are cooler in the fluid model, with a correspond- 
ing decrease in collision length and increase in (a e ). This 
reemphasizes the importance of developing a better model 
of the sheath boundary conditions for fluid simulations of 
the SOL. 

One might notice that the t = 100/is temperature pro- 
file from the non-local model in Figure 2 looks odd: it 
increases near the boundaries instead of decreasing mono- 
tonically from the centre of the domain. The reason is 
that as the density from the ELM reaches the boundary it 
causes a temporary glitch while the boundary conditions 
adjust, with a very steep velocity gradient appearing near 
the boundary. This causes the convective term in 10 to 
heat the electrons, resulting in the upturn in temperature. 
The diffusive model allows enough heat-flux to conduct 
away this extra heat, but the non-local model does not, 
which is why this behaviour is noticeable only in the non- 
local model. The reason for the glitch may be the way the 
ion viscosity is calculated. A viscosity limiter is needed 
since without it the viscous heating would increase the 
temperature of the SOL plasma above the temperature of 
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Figure 4: Spatially averaged value of the flux limiter, a ei 
that would match the value of the non-local heat-flux. 

the source, which is clearly unphysical. On the other hand 
without some viscosity the velocity of the plasma would 
increase to above the sound speed. This suggests perhaps 
that a better method of calculating the viscosity may be 
useful. It is possible to calculate the ion viscosity in an 
analogous manner to the calculation of the electron heat- 
flux used here[24[. The neglect of time derivatives used in 
this approach is obviously not strictly valid for ions when 
ELMs occur: the plasma does then evolve significantly 
on ion timescales. Nevertheless the calculation would be 
valid both before the ELM begins and in the steady state 
reached if the ELM continues long enough. It might there- 
fore be reasonably hoped that, while not quantitatively 
accurate, it would interpolate between the two states in a 
qualitatively better way than a limiter which can be set 
only for one or the other state. 

2.1 Model Convergence 

In the non-local model one can increase the number of 
moments used; the effect is to extend the validity of the 
model to longer collision lengths[16], though at the cost of 
computational speed. The simulations reported here used 
L = 20, K = 20 (400 moments). To verify that this num- 
ber is appropriate, we ran simulations with fewer (L = 10, 
K = 10; 100 moments) and more (L — 40, K — 40; 
1600 moments) moments. As shown in Figure 5, the 400 
and 1600 moment simulations agree well across the whole 
range of temperatures in this simulation; they do differ 
somewhat in the size of their response to the 'glitch' men- 
tioned above but since that is unphysical this difference 
is not relevant. At the highest temperatures (and hence 
longest collision lengths) the 100 moment version gives 
significantly different results. 

3 Implementation 

The computation of the non-local heat-flux requires many 
integrals to be calculated: one per moment per field line. 
Even though the integrand depends on the point of eval- 
uation, it is not necessary to compute separate integrals 
for each point. The dependence, equation (5), is simple 
enough that integral at each point can be calculated from 
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Figure 5: Comparison of the electron temperature profiles for 100 (solid, blue), 400 (dashed, red) and 1600 (dotted, 
green) moment heat-flux models: on the left at t = 15/us (as the 100 moment model begins to diverge) and on the 
right at t = 50/is (near the peak electron temperature). 



that at the point either to the left or to the right (de- 
pending on the sign of the eigenvalue associated with that 
integral) so that the integrals at every point can be com- 
puted in a single loop over the field line. The integrals 
for different moments and for different field lines are inde- 
pendent, allowing efficient parallelization by overlapping 
their computations on different processors to minimize idle 
time. Since this is a one-dimensional, field-parallel cal- 
culation, the scaling as the grid is divided in the other 
directions (i.e. radially and toroidally) is perfect (i.e. lin- 
ear) since processors dealing with separate field lines do 
not need to communicate. This favours the division of a 
three-dimensional simulation grid as finely as possible ra- 
dially and toroidally before starting to distribute the grid 
points on a single field line between different processors. 
However, other operations have diametrically opposed re- 
quirements: BOUT++ does not split the simulation grid 
at all in the toroidal direction in order to allow efficient 
computation of toroidal Fourier transforms; and perpen- 
dicular Laplacian inversions favour splitting of the grid 
parallel to the field rather than radially[19]. Therefore it 
is important that this heat-flux calculation scale efficiently 
when the grid is split in the parallel direction, but with 
several field lines on each processor. The sub- grid assigned 
to each processor will have at least as many points as the 
toroidal size of the grid, and likely several times that: this 
is the reason that Figure 6 shows the scaling for various 
numbers of field lines per processor, even though the most 
efficient splitting for this calculation would be to have just 
one field line per processor. The maximum number of field 
lines plotted is 64 since, for a grid with 256 parallel points, 
the evaluation time per grid point for larger numbers of 
field lines did not change significantly from that for 64. 
Figure 6 shows that scaling (for a 256 point grid) is near 
linear up to 16 = \/256 processors, falling to ~ n~ - 6 above 
that as the communication time becomes relatively more 
important. The efficiency (for 64 field lines per processor) 
at 64 processors, i.e. 4 parallel points per processor, is 56% 
compared to using a single processor. 

The grid for a typical three-dimensional simulation in 
BOUT++ might have 256 radial, 256 field-parallel and 
256 toroidal points; dividing the grid to put 4 radial and 4 



parallel points on each processor would mean using 4096 
processors. From the point of view just of this parallel 
heat-flux calculation this means 64 independent sets of 64 
processors, with each set having 1028 field lines. As the 
number of field lines is greater than 64, the efficiency will 
still be 56% and the total evaluation time (for the heat- 
flux calculation) for the 256 x 256 x 256 grid will then 
be ~ 2.4s. Including this non-local heat-flux in a three- 
dimensional simulation will not be detrimental to the scal- 
ing performance of BOUT++, and so we anticipate good 
efficiency up to several thousand processors, as has been 
previously demonstrated for BOUTH — h simulations [19]. 

The spatial convergence of the calculation compared to 
a 9192 point grid is shown in Figure 7. 256 point resolution 
has a root mean square relative error of 4.7 x 10~ 5 . The 
input data for the test are some functions chosen to have 
gradients comparable to the biggest ones found during the 
ELM simulations (modelled on the t = 2fi,s slice). 

4 Conclusions 

By using a non-local operator to calculate the electron 
heat-flux, we can account for some kinetic effects in fluid 
models, without needing to set parameters by comparison 
with kinetic simulations or experiment. This is especially 
beneficial when transients such as ELMs occur because 
the non-local heat-flux can respond self-consistently to the 
changing conditions in the plasma in a way that flux lim- 
iters cannot. 

In one dimensional simulations, the effect on machine- 
relevant parameters, such as the heat-flux at the divertor 
targets, is not significant as their behaviour is dominated 
by the ion dynamics (Figure 1) which, being largely con- 
vective, are well described by fluid models. 

However, there is a substantial change in the shape of 
the electron temperature profiles (Figures 2 and 3) which 
will alter the drive of turbulence in three-dimensional sim- 
ulations and thus will affect, for example, the width of 
the strike-point on the divertor. The principal advant- 
age of the technique described here is that it will allow 
self-consistent three-dimensional simulations: since the 
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Figure 6: Scaling performance of the non-local heat-flux 
computation for 400 moments on a 256 point grid on 
HPC-FF. Extra field lines allow more efficient overlap- 
ping of computations, improving the scaling compared 
to the single field line case. The black, solid line shows 
linear scaling and the red, dashed line shows n~ 6 scal- 
ing. 



Figure 7: Root-mean-square relative error compared to 
calculation with 9192 grid points. The red, solid line is 
oc -nT 1 and the green, dashed line is oc nT 2 . 



sources of heat and particles on each field line would then 
be determined by the cross-field transport, achieving a 
similar effect through flux limiters set by comparison to 
one-dimensional PIC simulations would be very challen- 
ging. 

The main limitation preventing a closer match of the 
electron dynamics between these non-local fluid simula- 
tions and PIC simulations such as those in [15] seems to 
be the response of the sheath-edge boundary conditions to 
kinetic effects during transients, which have not yet been 
included in our fluid simulations. 

The non-local heat-flux has been implemented in 
BOUT++ and parallelized efficiently. It is ready to be 
included in three-dimensional simulations of the SOL. 

4.1 Future work 

The next step in this work is to move on to three- 
dimensional turbulence simulations of the SOL, to invest- 
igate the effects on turbulence and perpendicular trans- 
port of the electron temperature gradients seen here that 
are missed by local heat-fluxes currently used in SOL fluid 
models. 

In order to correctly describe realistic tokamak geomet- 
ries, the model should be extended to include the effects 
of V.B on the heat-flux calculation. The effect will be to 
add an additional drive term to equation (1) which, as it 
depends only on the Maxwellian part of the distribution 
function, should add only a modest amount of additional 
complication. 

In the non-local approach one has more information 
about the electron distribution function than just the 
density and temperature. It may be possible to use this 
information to improve the modelling of the sheath and 
so find more accurate boundary conditions for SOL fluid 
models. Whether by this method or another, finding a way 
to include some kinetic corrections to the sheath trans- 
mission coefficients would greatly improve the accuracy of 



fluid descriptions of parallel transport in the SOL. 

Since convection dominates over conduction for the ions, 
kinetic corrections are less important; flux and viscosity 
limiters describe the behaviour of the ions better than is 
the case for the electrons. There may however be some 
scope for improvement in this area from application of a 
non-local model to the ion dynamics also. The variation 
of the ion viscosity limiter (inferred from a PIC simula- 
tion) through an ELM crash is not as dramatic as that of 
the electron heat-flux limiter[15]. However, the value of 
the viscosity limiter does have a significant effect on the 
dynamics, affecting in particular the speed reached by the 
plasma and hence the time for the ELM to reach the tar- 
get. Though a non-local calculation on the same lines as 
that used here for the electron heat-flux[24] is not guaran- 
teed to be valid for the ion dynamics due to the neglect of 
time derivatives, it seems worth investigating whether it 
might be able to give a better qualitative agreement with 
kinetic simulations and therefore be of use for improving 
fluid turbulence simulations. 
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